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We propose an improved fast multi-orbital impurity solver for the dynamical mean field theory 
(DMFT) based on equations of motion (EOM) of Green's functions and decoupling scheme. In this 
scheme the inter-orbital Coulomb interactions are treated fully self-consistently, involving the inter- 
orbital fluctuations. As an example of the derived multi-orbital impurity solver, the two-orbital 
Hubbard model is studied for various cases. Comparisons are made between numerical results 
obtained with our EOM scheme and those obtained with quantum Monte Carlo and numerical 
renormalization group methods. The comparison substantiates a good agreement, but also reveals a 
dissimilar behavior in the densities of states (DOS) which is caused by inter-site inter-orbital hopping 
effects and on-site inter-orbital fluctuation effects, thus corroborating the value of the EOM method 
for the study of multi-orbital strongly correlated systems. 
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I. INTRODUCTION 

The dynamical mean field theory (DMFT), which 
was proposed^ anc j developed in the past decade, is 
presently one of the most widely used methods for study- 
ing strongly correlated electron systems. Although it is 
exact only for infinite dimensional systems, it has been 
proven that it can provide a good approximation even for 
three dimensional systems. For low dimensional systems, 
a cellular DMFT has been developed in order to take into 
account spacial fluctuations. For a review of DMFT, see 
Refs.[3]and|U 

In the DMFT, the main tasks are to solve the so- 
called impurity problem and to derive the DMFT self- 
consistency conditions for a specified lattice of a real ma- 
terial. For the former task, several numerical techniques 
have been developed to construct the impurity solver, 
e.g., the quantum Monte Carlo (QMC) methodpthe ex- 
act diagonalization (ED) method,^ 7 numerical renormal- 
ization group (NRG) methodJ^ED density matrix renor- 
malization group (DMRG) method,^ the c ontinu ous 
time quantum Monte Carlo (CTQMC) methodpJEIl and 
the equation of motion (EOM) methodPHUl The EOM 
method, which is the focus of the present article, works 
directly on the real frequency axis and for any temper- 
ature, and in addition has numerical advantages as it 
is computationally faster than other numerical methods, 
because other methods are all computationally expensive 
due to the involved matrix manipulations, in which the 
Hilbert space dimension increases exponentially with the 
number of orbitals. Moreover, the EOM method is a very 
general method that can be implemented to a large vari- 
ety of systems in different cases. Also, through the study 
of equations of motions, we can learn more about the 
physics underlying interesting many-body systems, and 
do not only focus on the numerical technique and code 
development. The EOM method has actually been used 
by J. Hubbard in his corner stone workP^HIH anc i j^g corL _ 
tinually been studied by various people in studies of the 



single-impurity Anderson model (SIAM) and Hubbard 
model (see, e.g., Refs. 19 26). Here we study the impu- 
rity solver based on the EOM method within the DMFT 
framework. This study will also be beneficial for the 
development of further techniques, for example cellular 
DMFT for low-dimension systems beyond the standard 
DMFT, or development of techniques for directly solving 
the multi-orbital lattice model with tight binding. 

On the basis of equations of motion and decoupling 
techniques, an infinite U single orbital impurity solver 
has been developed in Ref. [T3] and a finite U single or- 
bital impurity solver with spin-orbital degeneracy N has 
been developed in Ref. 14. Recently, a multi-orbital 
equation-of-motion (MO-EOM) impurity solver has been 
developed by us in Ref. [T_3 which correctly captures the 
physics of many-body systems with the inclusion of the 
inter-site multi-orbital hopping interactions. However, 
in Ref. [T5]we have used the mean-field approximation for 
the on-site inter-orbital Coulomb interactions and have 
neglected the inter-orbital fluctuations, which is approxi- 
mately valid for those systems where the inter-orbital ef- 
fects are weak or there exists strong screening. For those 
systems with strong inter-orbital fluctuations, which can 
happen, e.g., when the orbital levels are energetically 
very close or when an applied external field exists, the 
inter-orbital fluctuations can be important and have to 
be included. 

In this paper, we have developed an improved EOM 
method for strongly correlated multi-orbital many elec- 
tron systems, in which we have implemented a higher- 
order decoupling scheme, where the on-site inter-orbital 
fluctuations associated with the inter-orbital Coulomb in- 
teractions are well included. A double check of the nu- 
merical results and iterative convergence procedure are 
introduced by calculating both the electron system and 
the associated 'hole' quasi-particle system in order to ob- 
tain accurate results. Using this MO-EOM scheme we 
find that the Kondo peak at the Fermi level is well re- 
produced for heavy fermions, and also the Mott metal- 
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insulator transition^ is observed. In addition, it is ob- 
served that the Hubbard bands are splitted due to the dif- 
ferent intra-orbital and inter-orbital Coulomb interaction 
strengths, and the quasiparticle DOS shows therefore an 
interesting multi-peak structure. This feature has been 
studied by tracing the difference between a two-orbital 
system with identical band widths for the two orbitals 
and that with different band widths. Comparisons are 
performed between results obtained with our MO-EOM 
method and those obtained with NRG and QMC meth- 
ods. 

The paper is organized as follows: The equations of 
motion, decoupling scheme, and physical assumptions are 
introduced in Sect.[n] Then in Sect. |III[ the results for the 
two-orbital Hubbard model as an example are discussed. 
Finally, a summary of the work is given in Sect. 
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II. DESCRIPTION OF MO-EOM METHOD 

We start our description from the model Hamiltonian. 
In many-body physics, the Hubbard model, in which only 
the hopping of electrons between the sites and the on-site 
Coulomb interactions are considered, is the simplest yet 
one of the most important lattice models. For a multi- 
orbital system, the Hamiltonian can be written as 

H = — ^ tijlm.f}i a fjmcj + ^ Uunufhill 
ijlma,i^j il 

~\~ ^ ^ UlnMj&'flilfjfl-irfKjf 
ilmacr' ,l<m 

i \^ (yi* ft f ,yl ft f \ (-1) 

~ / j V Ima J ima J ila ~ v Ima J ila J ima } t \ J 

ilma ,Z<m 

where i,j are the site indices, I, m are the orbital indices, 
a, a' are the spin indices, and f] l(T , fu a are the creation 
and annihilation operators for electrons with a spin in 
l-th orbital, respectively. The first term describes the 
hopping of electrons between the different sites. The sec- 
ond (third) term is the on-site intra-orbital (inter-orbital) 
Coulomb interaction term. The last two terms are the on- 
site inter-orbital single hopping of electrons, where V(* 
and V{ mcs are the inter-orbital hopping amplitudes for 
spin a between the l-th and m-th orbitals. 

In the dynamical mean field theory, the lattice model 
is mapped to an impurity model, usually the single- 
impurity Anderson model (SIAM), and the interactions 
between sites are mapped to the interactions between the 
impurity and a bath. For Eq. ([I]) , the mapped SIAM has 



the following Hamiltonian, 

n imp = S ^ C \ka C lka + ZflaflJla + ^ U llM^ 

kla la I 



Ika 

y , (Ylma fma fla + ^Ima fla fm< 



(2) 



lma,l<m 



The first (second) term is the energy of conduction elec- 
trons (localized electrons), where the electrons in differ- 
ent orbitals are labeled with the orbital index I. hi a = 
fiofia i s the occupation number for localized electrons 
with spin a in the l-th orbital, and Sfi a is the orbital 
level. The third term is the on-site intra-orbital Coulomb 
interaction term, and the fourth summation term is the 
on-site inter-orbital Coulomb interactions between elec- 
trons of the Z-th orbital and m-th orbital. The fifth sum- 
mation including two terms is the hybridization between 
the localized electrons and the baths. The last two terms 
are the on-site inter-orbital single hopping, where the site 
index has been dropped comparing to Eq. ([I]). 

In our studies, we have used the double time 
temperature-dependent retarded Green's function in 
Zubarcv notation,^ 



G AB {t,t') = <A(t);B(t') » 

= -ie(t-t')([A(t),B(t')} + ), 



(3) 



where A(t) and B{t') are the Heisenberg operators, and 
Q(t — t') is the Heavyside function. Applying the Fourier 
transform, we have obtained the Green's function in uj 
space, which should satisfy the equations of motion 



lo < A-B >= ([A,B}+)+ < [A,H imp ];B » 



(4) 



For a multi-orbital electron system, if we define the 
following anti-commutation relation for the operators, 

[fla,fma'} + = 0, [fl a ,f m(7 '] + = 0, 

[fla' fma'] + ^Im^aa' > [fla->f. mcr r] + ^Im^aa'i 

we can calculate with Eq. Q and obtain the first two 
equations of motion as follows (here we only show it for 
one orbital, m), 
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{bJ -\- jl & finer) < ^ I ma-, frna ^ > ^ ^tnrn < ^ ^mu' fma'i finer 

+ Yl ( V ™ k ° < C ™ fc -; t » - 1] < //ct ! /mo- » ), (5) 

(U) + fJ, - Ekma) < C mfc(T ; ffywj > = V mmfeo - < / mcr ; > + ^ Kn/fca < /i<7) /mcr >/#m, (6) 



where /i is the chemical potential. Moreover, we have em- 
ployed the notation a' ^ a which will also be employed 
in the following context. 

The first term of Eq. |5]) (i.e., 1) on the right hand side 
(RHS) reflects the existence of a particle with spin a in 
the m-th orbital itself. The second term on the RHS of 
Eq. ([5]) reflects the fluctuation of the spin a in the m-th 
orbital accompanied by the spin a' in the m-th orbital, 
which can be considered in another way as a fluctuation 
of the spin a when spin a' exists. The third (fourth) 
term gives the fluctuation of the spin a in the m-th or- 
bital accompanied by the spin a (a') in the l-th orbital, 
respectively. Furthermore, Eq. ([6| describes the situation 
that the electrons hop from the bath to the localized or- 
bital, in which event the electrons can hop to any of 
the orbitals of this impurity, while in Eq. ([5 ) there is the 
procedure that electrons hop from the loca ized orbitals 
to the bath. Therefore, we can imagine as a physical pic- 
ture that one electron hops from the m-th orbital of the 
impurity to the bath and then hops back to the m'-th 
orbital, where m! can be identical to or different from m. 
Consequently, the hybridization parameters V^ nmkrJ and 
V* mka label the processes that one electron comes from 
the m-th orbital and return to the m-th or i-th orbital, 
respectively. Due to the preservation of the charge, there 
exists the relation that V£ k(T = V* mmka + J2i 
Note that, because the SIAM is a mapped SIAM and 
therefore the bath in this mapped SIAM is actually a 
mapped virtue bath, it indeed corresponds to a real pro- 
cess in the lattice model in which one electron hops from 
one orbital on the studied site to any other orbital on any 
other site, which is the physical picture. One should note 
that V^ nmka and VL lfc(T will not appear in the impurity 
Hamiltonian Eq. pL see Fig. 1 in Ref. H5l 

From the derivation of Eqs. ^ and we can note 
that, when calculating the equation of motion of one 
Green's function <C A; B new Green's functions will 
appear, which are called higher order Green's functions 
than <C A;B ^>. If we calculate the equations of mo- 
tion of these higher order Green's functions, even higher 
order Green's functions will appear in these newly de- 
rived equations of motion and these equations of motion 
are also called higher order equations of motion having 
a higher order than those in the previous step. Each 
Green's function only has one order. Only those Green's 



functions that first appear and never appeared in previ- 
ous lower order equations of motion will receive its order 
by observing in which order equations of motion they 
first appear. Here the order of the Green's function ap- 
proximately labels the weight of the interaction associ- 
ated with this Green's function. Repeating this proce- 
dure again and again, more and more equations of mo- 
tion will be calculated and higher and even higher order 
Green's functions will appear, which is an infinite pro- 
cedure. Then a decoupling scheme will be employed to 
truncate this procedure and approximate those higher or- 
der Green's functions (higher than this truncation) with 
the product of the lower order Green's functions and rela- 
tive correlation functions. Thus the equations are closed 
and can be solved. The decoupling scheme gives under 
which order the interactions are treated exactly and all 
interactions higher than this order are treated approxi- 
mately. The higher this order is, the more accurate the 
scheme is. One can always get satisfactory precision of 
the results by choosing appropriate decoupling scheme, 
e.g., to choose one order higher decoupling scheme if the 
employed one is not sufficient. 



In Ref. [TSl we have treated the on-site inter-orbital 
Coulomb interactions with the mean-field approximation 
which may cause a loss of interesting information asso- 
ciated with the inter-orbital fluctuations. In the present 
treatment we take the inter-orbital fluctuations fully into 
account. We have previously noticed that, considering 
the pre-occupied charge besides the two spins in the 
two-particle Green's function that we are studying, the 
Coulomb interaction strength should be modified accord- 
ingly to be an effective onePS The DMFT will correct 
the shape of bands and the charge occupations, where 
the latter will cause a change of the Coulomb interac- 
tions between correlated electrons. Moreover, having the 
right Coulomb interactions will shift the bands to right 
positions, while the band positions and the shape of the 
bands will change the occupations. This implies that 
occupations, Coulomb interactions, band positions, and 
shape of bands have to be determined self-consistcntly. 
Now, since the inter-orbital fluctuation terms have here 
been taken into account, there should be three equations 
to describe the intra-orbital and inter-orbital Coulomb 
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interaction strength, respectively, 



where 



U mm = U mm + ^ (^mao-^ia + ^inW^ia') , (7) 



U, 



eff 
Imaa 



^ {Ul'-maa^Va + Ul'macr'fH' a') > 

I' 



(8) 



Ulmo-0-' ~\~ UfjimTlrfia-' -f~ Ul 7rL(7( jTli(j 
■ma a 

m< a + Ui> ma(T >ni /<7 /), (9) 



+ E w 



where I ^= m, I' ^= I and /' 7^ m. 

As a next step, we derive the equations in a form with 
the total hybridization functions A mCT , where for the sim- 
plicity we neglect the on-site inter-orbital direct single 
hoppings as we have done in Ref. [TS] and these on-site 
inter-orbital direct single hoppings will be treated in a 
forthcoming work. Therefore, we can obtain the follow- 
ing equation of motion, 



(UJ + fJb- Efmcr ~ A mCT ) X 

«/r„ CT ;/L»=i + ^, 



eJJ 



< n 



ma' fma 'i fma 



> 



+ E (Vim™ < nia-fma] f ma > 
+ U lmaa < flla'fma; fma > ), (10) 



A mm(T 



l,l^m 



LO + H — S m ka 



V m kaV lmka 



+ M - £ mfe CT < fma] fma > 



(11) 

(12) 
(13) 



Here the A mmcr (A; mo -) are the on-site indirect identical 
orbital (inter-orbital) hybridizations, which relate to the 
inter-site diagonal (off-diagonal) hopping terms. Under 
this assumption, there is the relation 



E V ™k<? < c mkt7 ; f ma >= A ma < f ma ; f ma > . (14) 



Now let us derive the EOMs for the higher-order 
Green's functions appearing in Eq. ^ and those newly 
generated equations, and then implement the decoupling 
scheme to those three-particle Green's functions to make 
the equations closed. For the sake of readability, we have 
put these EOMs in the Appendix. 

Solving the derived closed set of equations, we obtain 
finally the single-particle Green's function for spin a in 
the m-th orbital, 



fma] fma 



l + A 

) + B ( n la + lib) + C(ni a > + I U 



where 



w + (J, — £f mcT - A mCT — A(A m(T • I la + I 2a ) - J2i [ B[A. ma ■ I lb + I 2 b) + C (A ma ■ 7 lc + I 2c ) 



(15) 



^4 — U^J. 



M - £fma - U^'m - ^{Ui mcra ni a + Ulmaa'^la' + WuiaUlc 
I 

A m(7 A. maa A macr ^ l^ma' jAig' 4~ A[ a ) , 



b = u: 



eff 



imaa 



UJ + [I - Efma - UiJ na(y - (U mm n m a> + Vlma'a^lo' + Un ma n ma > + Un la rn a >) 

— ^^{Ul'maani'a + Ul'maa'^Va' + 2Uni> a ni> a') ~ A ma — A m ba — A m fc CT 
V 

-{nma' + nia')A ma - ^(n;v + fl Va i)A m 



(16) 



(17) 
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° - U lma' 



4 A 1 £frna ^Imcr'a (U m m n ma' 4 U[ m aa n la 4 Un mcr n mrT i 4 U ni a Tli a i^) 

— '^j^Vmaa^Va + 'Uvmaa'^Va' + 2t/n;' CT % CT /) — A„ lcr — A mccr — A mccr 
V 

-{n ma > 4 ni a )A ma - ^{n Vo 4 n ; , CT /)A m 



(18) 



_ Hmnfccr' Vmka' ^ y^ ^Imka' Vmka' *C fla' fma^ ^ Qg\ 

^ ''' + " + - p <"«x - £mb' „4^ W + H 4 £f m a> ~ £fma - £ m ka' < f m(T /; fj, > ' 



k 



W + /X + £/ roCT ' ~ £/ mCT - £ m fc<x' ;fe ^ m W + [L 4 £ fm „, - £ /mff - E mka > < / mff( ; /t^, > ' 

^ = y^ ^jjfcojjfcg |_ y^ ^ifcgjjfcg < /Wj /;t > ^0) 

-^-^ OJ 4- // 4- £fr- — Ef™. — £;i~, -^-^ W 4 (7, 4- E/u — E»™. — En. ^ • /t ^ ' 



A mC(J = V V Hk^Y^ h V" ^jW^W < /iv i /;t' > / 21 n 

t w + f + e /^' - e /w - j,^, wl/il e /J<r / - £/„ l(T - £; fc <x' < / J(T , ; > ' 

T/* V 

X \ v mka' v mmka' . 

^ ^ T /^T £mka' £ fma' £ fma U mma 

y^ Vmka'Vlmko' •C fla' j / mff ' ^ /^x 

fciT^n W + ^ + - e /™' - £ /™^ - ^n"™ < /ma' ! /L' > ' 



A m6(J = V VlkgVllka h y^ VjkgVl'lka < fl'a'Jla > ^ 

a _ \^ Viiv'Viik*' 

Y W + + £ "=:o-' - £fla> ~ £fma ~ Ummc 

y^ Vi* ka ,Vvika' < /;v ; /L > / 2 4) 

kVV^l UJ + t i + £ lk<y> - £fl*' - Sf-ma ~ U mrnc < fl(j' \ f\ a , > ' 



_ ^ Ymh^L (■fmcr' C, mka') Vmka' ( C mka' frag' ) ^ 

_ ^ _ Vmkcr'Vmka' ( c m k'cr' c rnka' ) Krifccr' ^mfca' ( c Lfcg-' c mk'a' ) -j 

CJ 4 4 £fmu' £fma £mka' CJ 4 /i 4 £mkcr' £ fma' £ fmcr Uj nma 



j ib _ y^ / Vlka(flcr C lk<j) Vlk<j{c\ ka fla) n 

j 2fc = y^ / _ ^L^fca( C L<T C 'fca) VlkaVlka{ c \ka C lk'a) \ ^ 

~f W 4 /l 4 £ fla - Sfrna - £lka W 4 /J 4 £;fc CT - £/ m(T - £ f la - U mrnb 

_ y^ /• ^ifcg' {Ila' c lka') Vmka' ( c lka' fla') ~j ^9) 

^ + £//o-' — £/ m<J — £lka' U) + H + £lka' — £ fla' — £ fma — U mmc 

j 2 _ y^ / _ Vlka'Vlka' (cjk'a'Clka') Vila' Vka' { c \ka' C ^'a' ) \ 

kk' w 4 /J 4 £fla' — £fma — £lka' W 4 jJ, 4 £lka' — £ fla' — £ fma — Ummc 
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where 

U mma = UgJ l + 2^(Ui nia + Ui ma o>'ni rT i + 2Um< T ni< 7 '), (31) 
i 

Uramb ~ ^Irnaa 2(£* ^ mmJ^rna' ~t~ ^ hmo '<7 '^la' ~t~ Un m( jJl ma ' ~t~ Ufli a Tli a i^) 

+2^2(Uv mcJa n v<y + Ui'maa-'fil'a' + 2Ufi V a 7lV a') , (32) 
/' 

Ummc — ^Imaa' nnri^ma' ~f~ ^Jlmaa^la ~t~ Uil rna Tl rrla f -)- XJ Tli a Tli(ji ) 

+2^(J7p mo . CT np (T + Ui'maa-mi'a' + 2Un Va n Va ,). (33) 
In the above equations, / 7^ Z'. Moreover, in the numerical calculations we will use the following Hermitian relations 



As we have mentioned in Ref. 1151 the hybridization 
functions in Eq. (151 will be obtained along with the 
DMFT self-consistency conditions. For the Bethe lat- 
tice, when one neglects the inter-site hoppings between 
different orbitals, the DMFT self-consistency condition 
should be 



A m(T t m <C fmcr't fma 



(34) 



When we take now into account the inter-site hoppings 
between different orbitals, we obtain that the DMFT self- 
consistency condition will be 



t 



L rn 
m ,2 



^ fma 5 fma ^ 



< fla'jL », 



(35) 



with t tot is the total amplitude summing over all the or- 
bitals, i.e., ttot — ^2 m tm- All the other hybridization 
functions can be obtained by interpolation. 

In this work, besides the inter-site inter-orbital hop- 
pings, we concentrate and fully include the inter-orbital 
fluctuations, because the Coulomb interactions are the 
most basic interactions and give the principal contribu- 
tions, and usually are also the only interactions incor- 
porated in most LDA+DMFT worksPHU other inter- 
actions, such as e.g. the spin-flip term and pair-hopping 
terrrP^, are considered to contribute less and contribute 
mainly to the shape of the Hubbard bands (for a detailed 
explanation, see Ref. [T5"|) . 



III. RESULTS AND DISCUSSIONS 

To investigate the accuracy of the derived method, we 
have studied the two-orbital Hubbard model in the half- 
filled, paramagnetic case with this MO-EOM impurity 
solver. The DOS for the quasiparticles are shown in 
Fig. [l] computed for the case where the band widths 




FIG. 1: (Color online) Quasiparticle densities of states for 
the halfRlled two-orbital Hubbard model on the Bethe lat- 
tice, computed with the MO-EOM method for the parame- 
ters: Di = D 2 = 1, T = 0.01, and J = 0. 



of the two orbitals are identical and the Hund's cou- 
pling constant is zero. It is observed that, along with 
the increase of the Coulomb interaction strength U, the 
system turns from the metallic state to the insulating 
state. The Mott metal-insulator transition occurs at 
nearly U = 1.45. We mention that the numerical cal- 
culation automatically fulfills the particle-hole symme- 
try and the sum rule that the integral over all the DOS 
equals identity. 

In Fig. [2] we show again the densities of states of 
the quasiparticles for the two-orbital Hubbard model 
at the halfhlling. Now we have chosen the parameters 
such that the two orbitals have identical band widths, 
i.e., Di = D 2 = 1, and the Hund's coupling constant 
J = C//4. The results for larger U reveal that the 
upper Hubbard band is in fact composed of three sub- 
peaks which are generated correspondingly by the intra- 
orbital and two types of inter-orbital Coulomb interac- 
tions. With the increase of U, these three sub-peaks are 
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FIG. 2: (Color online) Densities of states calculated for the 
halffilled two-orbital Hubbard model on the Bethe lattice for 
different U, with the parameters: D\ — D2 = 1, T = 0.01, 
and J = [7/4. 
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FIG. 4: (Color online) Quasiparticle densities of states for the 
halffilled two-orbital Hubbard model on the Bethe lattice with 
particle-hole symmetry at different U, with the parameters: 
Di = Z> 2 = 1, T = 0.01, and J = (7/4. 
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FIG. 3: (Color online) Computed densities of states for the 
halffilled two-orbital 'hole' system on the Bethe lattice for 
different Uhoie, obtained with the parameters: D\ — D2 = 1, 
T = 0.01, and J hole = U hole /4 = -[7/4. 



split more, due to the difference between the intra-orbital 
and inter-orbital Coulomb interaction strengths, i.e., the 
existence of the nonzero Hund's coupling constant. The 
result in Fig. [2] is obtained from the evolution with Eq. 
( 15 ), and it follows the sum rule. The word 'evolution' is 



here used for the iterative approach, iterating from an ini- 
tial guessed Green's function to a final converged result. 
Fig. [2] shows however an asymmetric DOS for the lower 
and upper Hubbard bands, something which is unusual 
and requires further investigations. As outlined in the 
following, this is related to how the iterative algorithm is 
performed. 

To study in more detail the issue of the particle-hole 
symmetry, we considering an electron system where all 
the states are initially fully occupied, the presence of a 
'hole' state (empty of electron) is equivalently a fermionic 



quasiparticle with spin index. Therefore, a many-electron 
system can be equivalently seen as a many-hole sys- 
tem, where the equivalent 'hole' system is similar to 
the electron system but only the interactions between 
the 'hole' states are with negative Coulomb interaction 
strengths. All our derived and obtained equations are 
also valid for these 'hole' states. Thus, we can imple- 
ment our equations of motion to this 'hole' system and 
set the hole position for <j spin channel in the m-th or- 
bital as -£ m fa- When J ho i e = U ho i e /4 = -Z7/4, we 
obtain consequently the DOS of the 'hole' quasiparti- 
cles, as shown in Fig. [3j However, the studied 'hole' 
system and the electron system are actually the same 
system, but are seen from different view angles. There- 
fore, the DOS obtained from the 'hole' system is also the 
DOS of the electrons. For any system, when the final 
self-consistently converged Green's function is obtained, 
no matter whether we calculate one next iteration as an 
electron system or as a 'hole' system, it will give identi- 
cal results, i.e., it will follow the particle-hole symmetry. 
One should note that this reasoning is not only applica- 
ble for the halffilled case, but for all arbitrary occupa- 
tions, i.e., all the systems should follow this rule. Com- 
bining this idea into the evolution, we have calculated 
both the electron system and the 'hole' system in our ge- 
netic algorithm schemeplwhich can also be considered to 
be a double-check procedure, and obtained the corrected 
DOS, shown in Fig. |4j for the halffilled two-orbital sys- 
tem with nonzero J = (7/4. The particle-hole symmetry 
is well recovered. We mentioned that both the lower and 
upper Hubbard bands show a splitting in sub-peaks for a 
large U due to the nonzero J. The lower Hubbard band 
is split, whereas in the DMFT self-consistency condition 
Eq. ( 35 ) for the Bethe lattice, the updated hybridization 



functions can not offer the information of this splitting. 
Combining the calculation of 'hole' states into the evolu- 
tion algorithm is a double-check process which can indeed 
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FIG. 5: (Color online) Quasiparticle densities of states com- 
puted for the halmlled two-orbital Hubbard model on the 
Bethe lattice at different J, for the parameters: D\ — D2 = 1, 
T = 0.01, and (7 = 4. 

safeguard the accuracy of the finally obtained result. In 
the genetic algorithm p2 this operation has the meaning 
to increase the diversity of the trial Green's functions, 
something which is important for the genetic algorithm 
so that the system will iteratively evolve in the right di- 
rection. If employed with linear mixing, the usual form 
of a searching scheme is 

G™ +1 = aGf w + (1 - a)G], (36) 

which means that in each iteration only a small amount 
of newly generated Green's function will be mixed with 
the Green's function of the last iteration so that the in- 
tegral equations are iteratively solved, a (and also the 
f3 below) is the mixing parameter and n is the iteration 
number. Our discourse above purports that now it should 
be modified to 

G n f +1 = aGf w + pGZZ + (1 - a - f3)G]. (37) 

In the genetic algorithm this does not cost additional 
time because in each iteration we will calculate a group 
of Green's functions. We can simply change some calcu- 
lations of electron system to 'hole' system. For details 
of the genetic algorithm for a continuous spectrum, see 
Ref. [3S] If the iterative search is performed with linear 
mixing, the required time will be doubled in serial cal- 
culation, but it is still much faster than other numerical 
methods. 

Fig. [5] shows the quasiparticle DOS at a fixed U but 
with different J, where the influence of the choice of J 
is studied. It can be observed that both the lower and 
upper Hubbard bands are split into sub-peaks by the 
Coulomb interactions and their splitting increases with 
the increase of the Hund's coupling constant J. Here 
we have chosen U = 4 which is a medium value. Taken 
the information in this figure together with that in Fig. 
[4j it is apparent that for a small U the Hubbard bands 
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FIG. 6: (Color online) Densities of states for the halffilled 
two-orbital Hubbard model on the Bethe lattice at different 
U, computed with the MO-EOM approach for two different 
band- widths, using the parameters: D2 = 2D\ = 2, T = 0.01, 
and J = 0. The DOS shown on the negative ordinate are those 
obtained for the wide orbital, those on the positive ordinate 
are obtained for the narrow orbital. 

can not be clearly split and will acquire only a change of 
the overall shape. For very large U, however, both the 
lower and upper Hubbard band can become split even 
into three peaks. 

In the above calculations, we studied the situation that 
the two orbitals have identical band-widths. Now let us 
turn to the case that the band-widths of the two orbitals 
are different, which can be occur, for example, for those 
systems that have both partial filled d and / electron 
shells, because the / electrons are considered to be more 
localized and should thus have a narrower band-width. 
To study this case we have used completely the same 
code, but only changed the parameters. Once a method 
can be generalized to more systems but without any ad- 
ditional work needed, it will be much easier to work with 
and will encounter more applications in LDA+DMFT 
calculations for real materials. In Fig. [6j we present the 
quasiparticle DOS for both the narrow and wide orbital, 
with the parameters such that the band-width of the 
wide orbital is twice that of the narrow orbital, and the 
Coulomb interaction strengths are identical, i.e., J = 0. 
For sake of clarity we plotted the DOS of the narrow or- 
bital on positive y axis and that of the wide orbital on the 
negative y axis. For the same U value, the narrow and the 
wide orbitals have different DOS. Also, we studied the in- 
fluence of U. With an increase of U, the narrow and wide 
orbitals simultaneously change both from being metallic 
states to insulating states. Hence, no orbital-selective 
Mott transition is observed. 

Fig. [7] shows the quasiparticle DOS for the param- 
agnetic two-orbital Hubbard model with different band- 
widths for the two orbitals, now computed with the pa- 
rameter J = (7/4. Similar features as shown in Fig.[6]can 
be observed, but at a higher critical value of U compared 
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FIG. 7: (Color online) Densities of states for the halfnlled 
two-orbital Hubbard model on the Bethe lattice at different 
U, computed with the parameters: D2 = 2D\ = 2,T = 0.01, 
and J = U /4. The DOS shown on the negative ordinate is 
that for the wide orbital, that on the positive ordinate for the 
narrow orbital. 
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FIG. 8: (Color online) Comparison of the densities of states 
for the halfhlled two-orbital Hubbard model on the Bethe lat- 
tice calculated with the MO-EOM to those computed with 
the NRG method for various U, using the parameters: D\ = 
D 2 = 1, T = 0.01, and J = [7/4. The DOS shown on the 
negative ordinate are those of the NRG method, those on the 
positive ordinate are computed with our EOM method. The 
NRG data are extracted from Ref. 36 



to the results in Fig. [6j which is due to the decrease of 
the effective Coulomb interaction strength caused by the 
nonzero J. In addition, the Hubbard bands display a 
splitting into sub-peaks in the large U region. 

Next, in order to address the accuracy of our MO- 
EOM impurity solver, we give some comparisons of our 
numerical results to those obtained with other methods. 
First, we present a comparison to the numerical renor- 
malization group method in Fig. [8] for the paramagnetic 
two-orbital system where the two orbitals have identi- 
cal band- widths, i.e., D\ = D2 = 1. The NRG data 



are extracted from Ref. [351 For the purpose of visibility, 
the NRG results are plotted on the negative y axis, while 
those of our MO-EOM method are plotted on the positive 
y axis. The quasiparticle DOS obtained with the same 
U are plotted in the same color. Overall, our MO-EOM 
method and the NRG method show a very good agree- 
ment no matter if it is on the widths and positions of the 
Hubbard bands or on the critical value of U for the Mott 
metal-insulator transition. In addition, our MO-EOM re- 
sult displays not only more micro-structures of the DOS 
caused by the nonzero J of the anistropic Coulomb inter- 
actions, but also shows a greatly reduced tail effect, which 
is an advantage of our EOM method that in our genetic 
algorithm evolution scheme the Lorentzian broadening 
can even be set as zero. With regard to these two points, 
our MO-EOM method has a computational advantage 
and is more accurate. 




FIG. 9: (Color online) Comparison of the densities of states 
obtained with the MO-EOM and QMC methods for the half- 
filled two-orbital Hubbard model on the Bethe lattice for dif- 
ferent U, with the parameters: D2 = 2Di = 2, T = 1/40, and 
J = (7/4. In panel (a), the quasiparticle DOS are compared 
between our MO-EOM method and the QMC method for the 
narrow orbital. In panel (b), the DOS are compared between 
our method and the QMC method for the wide orbital. In 
both panels, the DOS shown on the negative ordinate are 
from the QMC method and those on the positive ordinate are 
from our MO-EOM method. The QMC data are those from 
Ref.rjU 

Next, we provide a comparison of the MO-EOM quasi- 
particle DOS to values obtained with the QMC method 
in Fig. [9] for the case that J = £7/4. The QMC data are 
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obtained from Ref. [37J an d au parameters are taken to be 
the same. In panel (a), we compare the DOS of the nar- 
row orbital correspondingly obtained with our MO-EOM 
method and the QMC method, where the QMC data are 
plotted on the negative y axis and the MO-EOM data on 
the positive y axis. In panel (b), the DOS of the wide or- 
bital are compared between the QMC and our MO-EOM 
method, and similarly here the QMC and MO-EOM data 
are plotted on negative and positive y axis, respectively. 
To some extent, both methods have an agreement on 
the Hubbard bands, for the positions of the orbital lev- 
els. However, the two methods show a different critical 
value of U and different behavior when the system ap- 
proaches the metal-insulator transition, which has been 
found in a previous studyp^ too, and is here confirmed. 
The reason for the different behavior is that, in our MO- 
EOM method, we have taken into account the inter-site 
inter-orbital hopping effects (which are also mentioned in 
Ref. [32]), and the DOS obtained with the two methods are 
actually calculated with a different HamiltoniarP^^see 
Ref. [T5|) . The inter-site hopping reflects the interactions 
between the orbitals at the studied site and the orbitals 
of other sites. But, when all the sites are identical, due 
to the translation invariance, the inter-site inter-orbital 
hopping can equivalently be seen as an interaction be- 
tween one orbital and all other on-site orbitals, whose 
influence will suppress the orbital-selective Mott transi- 
tion. Moreover, we have taken into account the inter- 
orbital fluctuations. For spin a in the m-th orbital, the 
spins existing in other orbitals have a similar contribution 
as spin a' in the m-th orbital, which can be seen from 
our equations. In Fig. [2] we have noted that, for the two- 
orbital paramagnetic system, the Hubbard band for one 
spin is actually composed of three parts corresponding 
to the other three spins. Accordingly, the Kondo peak 
should also be composed of three parts. This means that, 
once there is a Kondo peak for one spin, there will be a 
Kondo peak for another spin. Therefore, the existence 
of the on-site inter-orbital fluctuations will suppress the 
orbital-selective Mott transition, too. We also note that 
in our MO-EOM result the Kondo peak is nearly pinned 
at the Fermi level, which is an important aspect for an 
impurity solver because in theory the Kondo peak should 
be pinned at one point at the Fermi levelpS while the 
QMC method does not show this feature. 



IV. SUMMARY 

In this paper we have proposed an improved fast 
and applicable multi-orbital impurity solver for the 
DMFT based on the equations of motion method with 
fully including the inter-orbital fluctuations beyond the 
mean-field approximation treatment to the inter-orbital 
Coulomb interactions. Our derivation has provided an 
alternative, fast and more precise way to obtain the local 



single-particle Green's functions. 

The derived impurity solver works directly on the real 
frequencies axis and for any temperatures. Comparing 
with other methods, it is less memory and cpu time ex- 
pensive. It has also given improved detail in the micro- 
structures of the DOS and shows an improved perfor- 
mance in reducing the tail effect. Furthermore, it still 
keeps the generality of the parameters. From our investi- 
gations it appears that our MO-EOM method can be gen- 
erally applied to a large range of multi-orbital systems, 
no matter if these have identical or different band-widths, 
nonzero J or not. In the present work we only studied the 
paramagnetic case, but the MO-EOM method can also 
be applied for the magnetic case. The comparisons of 
our numerical results with those obtained with NRG and 
QMC methods attest that our method shows not only a 
good agreement but also more precise micro-structures 
of the Hubbard bands, and also some distinct differ- 
ent behavior in approaching the Mott metal-insulator 
transition with the self-consistent inclusion of the inter- 
site inter-orbital hoppings and the on-site inter-orbital 
Coulomb interactions. Therefore, our MO-EOM method 
appears to be a good candidate for a fast impurity solver 
for the dynamical mean field theory and for LDA+DMFT 
calculations of multi-orbital strongly correlated electron 



systems 
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Appendix A: Equations of motion for the 
two-particle Green's functions 

In order to make the paper more easily readable, we 
present here the equations of motion for the higher-order 
two-particle Green's functions that are needed to derive 
the final single-particle Green's function Eq. ( 15 ) given 
in this paper. In addition we also list the decoupling 
scheme for the three-particle Green's functions appearing 
in those equations of motion. 

The EOMs of the higher-order Green's functions are 
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(w + n - e fma )G™ f = n m<J < + U mm < h ma ,f mi7 \ fl ia > 

ma' < ni a >n 

ma' f men fla ) 

i 

^mfefr' ^ ^mka' fma 1 fma , /m<r ~~^Vrnka ^ Tlma' 'C-mka ) /m<r 

+ V;„ fcCT ' < fma' c mka'fma\fma > ) ( A1 ) 
(W + /i - £ /mCT )G^™ / = nj ff + £/jWa < fllafma; fma > 

+ < hma>hi a f ma ; p ma > + £W< <J < h ltT ,hi a f ma ; p ma > ) 

+ 51 (E^Wcr < ni a niaf ma \f ma > i^'maV < n v a ,fii a j ma ; f} na > ) 

r 

+ X] ( ~~ ^ <C C lkaf^fma; fLa > +Knfea < ?WC TO fc CT ; fla > 

+ VJ fcCT <// ff c tt(7 / TO(7 ;/^ CT >) (A2) 

(CJ + /i — £/ mcr )G^/y = n; CT ' + Ulma'u <C hla' fma', fma ^ > 

~t~ {Umm < hma'hl a > f ma ; flna > + £ / im<x C r < fll a h W f m a', fjna > ) 

+ X (U Vmaa < haha'f m a\fma^+ U l'rna'a^h'a'nia'f m a\f m a^ > ) 



V 



Vlka' ^ C jfc<T' /if' /m<r i /mcr + Knfcer ^ ^ma' c mka', fma 



/ , V Hfecr' ^ Ika 
k 



+ Vlka' « ft,Clka, f mrT ; f ma » ) (A3) 

(W + /! - £ mfccr )G™ = V^ TOfecr < n mCT '/m<r; ffna > + X (^mkV < fla' C mk' a> Cmka', fma > 

fc' 

^ Knfe'cr' <C cj^j,,^, fma'Cmka', fma ^ ) + X/ ^Imka ^ ^ma' fla', fma ^* (A4) 



(W + Ai - e mfeCT )G^, c = Cmfe<r < hlv'fma] f ma > + ^ (^'t' < fl> c lk>a>C m ka', fma > 

fe' 

- Vjfc'<x' < J lk , a ,fla'Cmka\fma > ) + X ^*mfea < /*'<x ! > (A6) 

(cj + /i + Efmcr' — e m ka' — s fma)Gfcf = (f mcr /C m ka') + V m mka' ^ T^ma' fma', fma ^ 

+ ^ ] ^Ymk'a ^ f mcr i c mka'Cmk'a', fma ~ ^rnfe'er' ^ C mfc'<r' c mka' fma ', fma ) 
fc' 

+ X ^mfca' « fla' fl-' fma', fla ^ (A7) 
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k' 

- Vj feV < 4, CT c ;feCT / ro<T ;/L>) + ^ V?, fc(T « /il/i'Jm,; /L > (A8) 

Z',2'^i 



(cj + /i + e /(CT , - e ;feff , - e fma )G l p cf = {f} a ,ci ka ,) + V^ felT , < h la ,f mr7 ; f ma > + ^ (V^ fc , ff < f mCT ,c mk(T >c mk > a ; fLo > 

fe' 

»)+ E ^ « /U'"/™ 4. » (A9) 



(cj + /i + £ m fe CT ' £f ma i — £f ma )G™ff — (clnka-i fma') + U mm c lnka' fmcr' fmcr', fmcr ^ 

C mka' fmcr' fma , /mo- ~^rUl m a'a *S ^Zcr' c m ka' finer' fma, fma ) 

/ 

— Knmfcff' <C flma'fma', fma ~l~ ^ ] (^mfc'er' C mka' C mk'a' fma', f m a S> 

fe' 

+ Knfc'rr' < cLka'fma'Cmk'cr-jLa » ) ~ E < f la' fma' fma', f La ( A1 °) 

l,l^m 



{uj + H + £l ka -Efl a - e f mcr)G l c y af = (c] k Jla) + U lmaa C^Jlafmcr-jLcr^ 

( ^ i" i" ^ i" "I - *\ 

+ (Umm ^ ^mcr' C-i k <jflcr fmcr] fmcr "^UmXa'a ^ flma' ^l ku fla fmcr, fmcr ) 
+ (Z7 H < hl<j>c] k(y f la .fmc,; fLcr > +£4n<7'<7 < flla'4 ka flafma', fmcr > ) 

'merer <C n.J' ac\ kcr flcr fmcr', fmcr +t^i'm<j'cr <C filer' c[ ka flcr fmcr', f mC r 



+ £/*'W < fli'crc\ ka flcr fmcr',. fmcr > +^'Zo-' C r < "^■ler'^lkcrfler fmcr', fmcr ) 

C \kcr C lk' a fmcr', fmcr 

fe' 

+ CfcV< 4 k J la Cmk'cr',fL^>)- E ^'^« f\'Jlorfmcr',fmcr~>, (AH) 



(W + M + EifefT' - £/Zcr' - £fmcr)G l ^ a , f - (4 klT , flc') + Ulmcr'cr < c\kcr'fl<?' fmcr', flier > 

+ (Umm <C ^mcr'C^kiji flcr' fmcr', fmcr ^Umlcr'cr' <C flma' c lka' fler' fmcr', fmcr ) 

+ (t/jj < nicrcjkcr' fl°' fmcr', fmcr > + £ / im<x<T < hlcrc\ ka ,flcr' fmcr', fma > ) 

+ ^ ] (Ul'mcrcr *C ttl' cr c \ k(J i flcr' fmcr', fmcr ^> +Ul'mcr'er *C ^ler' c } klJ i flcr' fmcr', f m a 

+ t^i'icro-' < ni'crc\ k(T ,fla'fmcr; fLe, > < filer' c] k<J , ,flcr> .fmcr', fmcr > ) 

- V HfeCT ' < filer' fmcr', f ma > + E < 4ka' C lk' a' fma', f m a > 

fe' 

+ CfeV< 4a'/^' C ™feV;/L>) - E fia'fl°'fma',fia^>, (A12) 

where the Coulomb interaction strength will change accordingly for the multi-orbital case due to the pre-existence of 
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charge, and we have used the following abbreviations for the Green's functions on the left hand side: 

^nf ^ f^mcr' fma, fma ^~*fcf ^ "ma' C " m ^a l fma ! fma 

@nc =< ^ fhntr' ' c mkal fma -^i ^~*cff =< ^~ C mka' f m a' fma, /m<7 • 

^naf =< ^ ^lafmat fma ^) ^facf fl a c lka fma', fma ) 

^nac~^ f^laCmka't fma ^> ^cfaf =< ^~ C lkafl< J f' mcr ''fma^ > - 

^na'f = ^ ^la' fma', fma ^Ja'cf =< ^~ fla' C lka' fma', fma ^ , 

G n a'c =< ^ nia'Cmka', fma ^cfa'f =< ^ C lka' fl<r' fma") fma ■ 

In the above equations of motion, we will drop the Green's functions involving odd number of operators of one 
arbitrary orbital, e.g. <C f}>a'fia' fma', f m a ^ etc., and attribute them to the effects of inter-orbital hybridizations 
which will be taken into account through the DMFT self-consistency conditions. 
For the three-particle Green's functions we used the following decouplings, 

< ni 

a^ma' fma] fma -^ > ~ n la *^ ^ma' fma \ fma -^ > ; (-^13) 

^ma' fma j fma -^ > ^ ^Zcr 7 < ^ ^ma' fma-, fma (■^^^) 

T^ma'Tlla fma\ fma --^ > ~ ^mcr' *^ ^la fma'i fma -^ > ? (-^15) 

< nia'tilafmcrlfha >~ < "iff /ma? /L ( Al6 ) 
Tlrna''N'la' fma i fma -^ > ^ T^ma' < nia>fma',fia^>, (A17) 

< haha' fma) fma >~ < hla>fma\f ma > • ( Al 8) 
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